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ABSTRACT 

We present a self-consistent model of a protoplanetary disk: ANDES' (AccretioN disk with Dust Evolu- 
tion and Sedimentation'). ANDES is based on a flexible and extendable modular structure that includes 1) a 
1+1D frequency-dependent continuum radiative transfer module, 2) a module to calculate the chemical evo- 
lution using an extended gas-grain network with UV/X-ray-driven processes surface reactions, 3) a module 
to calculate the gas thermal energy balance, and 4) a 1+1D module that simulates dust grain evolution. For 
the first time, grain evolution and time-dependent molecular chemistry are included in a protoplanetary disk 
model. We find that grain growth and sedimentation of large grains to the disk midplane lead to a dust-depleted 
atmosphere. Consequently, dust and gas temperatures become higher in the inner disk (R < 50 AU) and lower 
in the outer disk {R > 50 AU), in comparison with the disk model with pristine dust. The response of disk 
chemical structure to the dust growth and sedimentation is twofold. First, due to higher transparency a partly 
UV-shielded molecular layer is shifted closer to the dense midplane. Second, the presence of big grains in the 
disk midplane delays the freeze-out of volatile gas-phase species such as CO there, while in adjacent upper 
layers the depletion is still effective. Molecular concentrations and thus column densities of many species are 
enhanced in the disk model with dust evolution, e.g., C0 2 , NH 2 CN, HNO, H 2 0, HCOOH, HCN, CO. We also 
show that time-dependent chemistry is important for a proper description of gas thermal balance. 
Keywords: accretion, accretion disks - circumstellar matter - stars: formation - stars: pre-main-sequence, 
astrochemistry 
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1. INTRODUCTION 

The planet formation and, in particular, the origin of the So- 
lar System are among the most fascinating astrophysical prob- 
lems that are far from being fully understood. The quickly 
growing number of detected exoplanets hints to ubiquitous 
planet formation in our Galaxy. Space-born facilities (e.g., 
Hubble, Spitzer, Herschel) as well as ground-based observa- 
tories (e.g., VLT, Keck, Subaru, PdBI, IRAM 30-m, SMA, 
early ALMA) provide unique information on the appearance, 
structure, chemical comp osition, and evoluti o n of nearby pro- 
toplanetary disks (e.g.. lDutrevetal.il 19971 iFukagawa et al.l 
20041: lAndrews & WiUiamsl 120051: (Hernandez et al.1 | 2007[ 
Natta et al l 1 20071: ISemenov et all l2010t ISturm et al.l bold 
Muto et alJ 12012b . Relatively compact sizes of ~ 100 — 
1000 AU and low masses of ~ 0.01 M Q make disks a chal- 
lenging target for observational studies. 

Another obstacle to investigate the formation of planets is 
an enormous range of physical conditions encountered in a 
protoplanet ary disk and a wide va riety of interrelated pro- 
cesses (e.g., William s & Ciezal201 ll) . The combined action of 
these processes defines the appearance of the disk in scattered 
light, dust continuum, and atomic and molecular lines. Mod- 
eling of continuum and line radiation implies knowing stellar 
spectrum, dust density, dust temperature, and size distribution 
as well as gas density, gas temperature, and molecular content 
throughout the disk, and in full 3D. If all this information is 
available, a multi-dimensional radiation transfer (RT) model 
can be used to build a synthetic disk map at any wavelength 
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(e.g., IWhitnev & Hartmannlll992t IMen'shchikov & Henningl 
ll997tlWolf et alJll999tlDullemond et alJl2002h . Due to com- 
putational difficulties to follow global disk evolution in 3D- 
MHD, particularly, coupled with chemical kinetics models, 
and the lack of necessary constraints related to the magnetic 
field structure, turbulence, grain size distribution, etc. , a disk 
model needs to be simplified. One can steadily approach the 
warranted level of physical complexity by adding new com- 
ponents to the model (e.g., going from ID to 2D geometry or 
from gray to non-gray radiative transfer) and comparing with 
observations at each development step. 

A num ber of disk models has bee n developed over time (see 
review in Dulle mond et a l. 2007b). These models have been 
based on an RT-based disk structure (either ID, 1+1D/2D, 
or 3D), molecular abundances, and dust and gas thermal 
balance. Disk models with detailed vertical structure and 
thermal balance regulated solely by dust heating and cooling, 
and, in so me cases, ac cretio n heating, have been d e velope d 
bv. e.g.. iBelletalJ (fl997h: IChiang & Goldreichl (fl997l): 
IMen'shchikov & Henningl (119971) : iDullemond & Dominikl 



(2004): iHueso & GuilloT (120051) . It has been typically 
assumed in such studies that the dust is well mixed with the 
gas, and its properties do not differ from properties of the 
ISM dust. One of the most widely used models of t his kind 
has been developed bv lD'Alessio et all d!998l 119991) . It has 
been extensively used in many subsequent studies as a tem- 
plate of the disk d ensity and temperatu re distributio n (e.g ., 
IChiang et al.ll200lt ISemenov et alJl2004t iFurlan et al.l 2006). 
Other similar models, utilizing more accurate frequency- 
dependent RT algorithms or other improvements (e.g., a 
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full 2D geometry, evolving disk stru cture, more r ealisti c 
dust opaciti e s) hav e been presented by Mal bet et al.l d2001l) ; 
iDullemondl d2002h : iNomural J2002h ; iGorti et al.l (120091) . to 
name a few. 

An important development of the protoplanetary disk mod- 
els was to account for the energy balance of dust and gas 
separately in dilute disk regions. There the rate of gas- 
dust collisions drops so low that the gas becomes ther- 
mally decoupled from the dust (e.g.. iJonkheid .et al.l l2004t 
iKamp & DullemondN2004l : IGorti & HollenbacHl200ll . The 
most recent and most advanced ad dition to this family , the 
"ProDiMo" model, is presen ted bv IWoitke et al.l (120091) and 
updated in iThi et al.l d2011l) and lAresu et all (1201 ll) . This 
model is based on iterative calculations of a 1+1D verti- 
cal hydrostatic disk structure, 2D frequency-dependent dust 
continuum RT, gas-grain and FUV-photochemistry to cal- 
culate abundances of molecular coolants, and an escape 
probability method to model non-LTE heating and cool- 
in g of the gas. It is derive d from fhermo-chemical model s 
oflKamp & Be rtoldil fcOOOh. IKamp & van Zadelhof fi 120011) . 
and IKamp & Dullemond ( 12004) . Since 2011 it includes 
X-ray-driven chemis try and heating via H2 ionization and 
Coulomb heating (lAresu et al.l 1201 1|) . Uniform dust abun- 
dances and power- law size distributions are typically assumed 
( lAresu et al.l 12012b , with opacities for a dust mixtur e calcu- 
lated by Effective Medium Theory (Bruggem anl 19351) . Abun- 
dances of molecules are calculated assuming chemical equi- 
librium and element conservation, which m ay not be a valid 
appro ach to disk chemi cal evolution (e .g., Barshav & Lewisl 
1 9761 lllgner et al.ll2004t ISemenov & WiebeCOl ll) ! 

Recent observations at IR and mm-/cm-wavelengths have 
shown that many disks around young stars of ages > 1 Myr 
have already a deficit of of small grains in the inner regions, 
r < 10-50 AU and the presence of large, pebble-sized dust 
grains in the midplanes compared to the pr istine ISM dust 
(e.g., Williams & Cieza 201 1; Williams 2012). From the anal- 
ysis of SEDs at millimeter and centimeter wavelengths, grain 
sizes of at least 1 cm have been inferred for many disks (e.g., 
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Rodmann et al. 2006; Lommen et al. 
2010t iMelis et alJl201 U iPerez et all 
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(1201 ll) have used high-resolution interferometric PdBI obser- 
vations to discern dust emissivity slopes at millimeter wave- 
lengths in a sample of young stars. Their analysis has shown 
that in the Taurus-Auriga star-forming region some disks 
show very low dust emissivity indices in the inner regions, 
characteristic of grains with sizes of > 1 mm, and slopes 
that are indicative of smaller grains toward the disk edges. 
In addition, Spitzer IR spectroscopy of silicate bands at 10 
and 20 /im has revealed efficient crystallization and growth of 
the sub-micron-sized ISM grains in warm disk atmospheres in 
many young systems, regardless of their ages, accretion rates, 
and d i sk masses (e.gJKess ler-Sila cci et al.l20 06: Furl a n et al] 
20091: Uuhasz et alJl2010t iMcClure et al.l l2010t lOli veira et al.1 
201 U lSicilia-Ag uilar et al.l201 ll) . The dust settling associated 



with grain growth reduces disk scale heights and flaring an- 
gles, and thus leads to less intense mid-IR disk emission than 
expected from conventional hydrostatic models with uniform 
dust, in accordance with observations of most T Tauri stars 
( Willi ams & Ciezall20Tl . 

As dust is a very important ingredient of the disk physics, 
evolution of its properties should also be considered in disk 
models. Usually both the grain growth and sedimentation are 
accounted for in disk models in a parameterized way, by as- 
suming an increased upper limit of grain size a max and arti- 



ficially changing the dust density and the slope of dust size 
distribution in variou s disk regions. For exa mple, expanding 
on their earlier works, D' Ales sio et al.l (12001 ) have studied the 
influence of dust evolution on the disk structure and its spec- 
tral energy distribution (SED). Grain growth has been simu- 
lated as an increase of a max up to 10 cm and change of the dust 
size distribution slope p from -3.5 to -2.5. In these models 
dust has been assumed to be well-mixe d with the gas. 

To study the effect of dust settling, ID' Alessio et al.1 (120061) 
have included two dust p opulations in the model , with differ- 
ent spatial distributions. D'Alessio et al. (2006) shown that 
the evolved dust model better reproduces observed millime- 
ter fluxes and spectral slopes. A similar approach to study 
the effect of dust settling on th e disk thermal and ch emi- 
cal structure has been taken by Jonk heid et al.l (120041) and 
iFogel et all (1201 ll) . Settling has been simulated using vari- 
abl e dust/gas mass ratio. A va riable a max value has been used 
by Aikawa & Nomura (2006) to investigate changes in disk 
density, gas and dust temperature, and molecular abundances 
due to dust growth. 

More accurate methods to model dust growth are mainly 
based on solving the coagulation (Smoluchowski) equation. 
Here the main attribute of the model is whether the dust 
evolution is computed for a fixed disk structure or the dust 
evolution and disk structure are mutually consisten t . Th e 
first approach is used , e.g., i n iNom ura & Nakagawa (2006); 
Schra pler & Henningl (120041) ; iTanaka et al.l (120051) : ICieslal 
(120071) , who used parame terized disk structure . The second 
approa ch has been use d bv | Schmitt et al.1 (ll997l):lTanaka et al.1 
(2005); N omura et al.1 (120071) : iTann irkulam et al.1 (l2007l) ~ 

An efficient scheme to tackle the modeling of dust coagula- 
tion, fragmentation, sedimentation, turbulent stirring around 
a 'snow line ' in a p rotoplanetary disk has been proposed by 
iBrauer et al.1 (12008). They have found that major factors af- 
fecting grain evolution are trapping of dust particles in gas 
pressure maxima and the presence of a turbulently qu i escent 
'dead zone' in disk inner midplane. iBirnstiel et al.l (120101) 
have updated this model by considering time-dependent vis- 
cous evolution of a gas disk. They have found that dust prop- 
erties, gas pressure gradients, and the strength of turbulence 
are more important factors for dust evolution than the initial 
conditions and the early forma tion phase of the protoplane- 
tary disk. Birn stiel et al.l (1201 ll) have shown that, upon evolu- 
tion, grain size distribution reaches a quasi-steady state, which 
however, does not follow the standard MRN-like power-law 
size distribution and is sensitive to the gas surface density, 
amount of turbulence, and disk thermal structure. 

The next step in proto planetary disk modeling was made 
bv lVasvunin et al.l d201 ll) . where detailed dust evolution was 
considered along with comprehensive set of gas-phase and 
surface chemical reactions. However, to calculate disk ther- 
mal structure, they take into account only two heating sources, 
namely, viscous dissipation and dust grain irradiation by 
the central star. It was shown that column densities of 
some molecules (like C2H, HC2„+iN (n = 0-3), H2O and 
C2H2/HCN abundance ratio) can be used as observational 
tracers of early stages of the grain evolution in protoplanetary 
disks. 

In this paper, for the first time, we consider the influence 
of dust evolution on the disk structure by combining the de- 
tailed computation of the radiation field with the dust growth, 
fragmentation, and sedimentation model. When computing 
the disk density and temperature we take into account the full 
grain size distribution as a function of location in the disk. 
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Gas temperature and dust temperature are computed sepa- 
rately, with taking into account the disk chemical structure. 
These t wo factors represent a m ajor improvement in compari- 
son with lVasyunin et all (1201 l)'s model. Also, a new detailed 
RT treatment is implemented with high frequency resolution 
from ultraviolet to far infrared. The organization of the pa- 
per is the following. In Section 2 the disk model 'ANDES" 
(AccretioN Disk with Dust Evolution and Sedimentation) is 
described. In Section 3 we present a physical structure for 
a typical protoplanetary disk computed both for both pris- 
tine dust and for evolved dust. Also, the chemical structure 
is described in this section, and specific features of the disk 
chemical compositions are presented for various dust mod- 
els. Discussion and conclusions follow. Details of gas energy 
balance processes and benchmarking results are presented in 
AppendixlAlandlB"] 

2. DISK PHYSICAL STRUCTURE 

A multitude of processes (gas dynamics, dust evolution, en- 
ergy transport processes, chemistry, etc.) makes modeling of 
protoplanetary disks a challenge. With the current level of 
computing resources a global 3D radiative MHD simulation, 
including gas and dust evolution and ch emical kinet i cs, re- 
mains a topic for the future (but see, e.g., iFlock et ail (12012b 
for such models). Nevertheless a sufficient understanding of 
protoplanetary disk physics may be achieved by detailed mod- 
eling of primary processes that govern its structure and obser- 
vational characteristics, and simplified modeling of secondary 
process es. This makes the problem t ractable. For Class II 
objects (Lada 1987; Evans et al. 2009) it is usually assumed 
that a disk structure is in a steady-state regime over a time 
span of ~ 1 Myr. This is supported by observations of disk 
kinematics via molecular lines and disk surface densities via 
(sub-)millimeter dust emissivity. The line profi les are indica- 
tive of Keplerian motion in most of the disks dKoerner et alj 
[1991 iGuilioteau & Dutrevll 19981: iPiitu et alj|2007l) . The es- 
timates of the disk masses and density distribu tions show 
that s elf-gravity is negligible for Class II objects (Isel la et ail 
2009). The assumption that these disks evolve on a diffu- 
sion timescale and not on a hydrodynamical one allows setting 
aside hydrodynamical simulations and reducing a 3D problem 
to a 1+1D problem. The azimuthal dimension is eliminated 
due to the axial symmetry of an unperturbed disk. The other 
two dimensions are usually split into the radial structure that 
is determined by diffusive evolution, and the vertical struc- 
ture that is derived from the hydrostatic equilibrium equation 
(iD'Alessio et al.ll998fclDullemond et alJ2002h . The 1+1D de- 
scription is suitable for dust continuum radiation transfer. For 
disk regions outward of a few AU a radial optical depth for a 
location close to the midplane is higher than the vertical op- 
tical depth, so that the dust temperature is mostly determined 
by vertical diffusion of radiation. A gain in computation time 
that is acquired by a ID radiation transfer, compared to a 2D 
RT, allows better frequency re solution, which is important for 
dust temperature calculations (Dullemond et al. 2002) and for 
modeling photochemistry. 

In this paper we adopt a 1+1D approach to calculate disk 
density and temperature. As the disk consist of two main in- 
gredients (dust and gas), the overall problem is reduced to 
calculating four physical quantities: dust and gas tempera- 
tures (7^,7^), and dust and gas densities (pd,Pg)- This allows 
to split the disk model into four blocks, calculating the corre- 
sponding quantities at each disk radius R: 



I. Dust temperature. Dust temperature and radiation 
field /„ are found by solving the radiation transfer 
problem in vertical direction. The following quanti- 
ties are considered as input: the dust density, its op- 
tical properties (absorption and scattering coefficients 
K v ,a u ), external irradiation and all necessary parame- 
ters describing non-radiative dust heating functions. 

II. Gas temperature. To determine the gas temperature, 
we solve the local energy balance equation, account- 
ing for various heating and cooling processes. Since 
gas heating and cooling rates depend on abundances 
of main heating/cooling species and their level popula- 
tions, it is necessary to include chemical reactions and 
simplified line radiation transport in the gas tempera- 
ture calculation. 

III. Dust density. The dust evolution is an essential part of 
our disk model. The surface density of dust is assumed 
to be equal to 1% of the total gas density, whereas 
its detailed vertical structure and size distribution are 
determined from the dust growth and sedimentation 
physics. We consider coagulation and fragmentation of 
dust grains and their redistribution due to turbulent stir- 
ring and gravitational settling to the midplane. We also 
consider disk structure with for comparison. 

IV. Gas density. We assume that the gas vertical structure 
is defined by the local hydrostatic equilibrium. In this 
case the gas density p g can be found if its temperature 
T g , mean molecular weight \i, and surface density E s 
are known. The surface density is assumed to be given 
by the predefined function £g(R). 

As all these quantities are not independent, we iterate be- 
tween the modules until convergence is reached. The overall 
computational flowchart for ANDES is shown in Figure Q] 
Assuming the surface density profile, we calculate dust evo- 
lution for 2 Myr starting from the MRN initial distribution. 
The resultant dust structure is then used to derive radiation 
field and gas disk structure using radiation transfer, energy 
balance and hydrostatic structure modules. As a fiducial dust 
model we also consider pristine grains with the following pa- 
rameters: 0.1 /jm in size, astronomical silicate, dust to gas ra- 
tio 0.01 at every location in the disk. The list of basic as- 
sumptions is: (i) a disk is quasi-static, axially symmetric and 
treated in 1+1D approach; (ii) the gas surface density is as- 
sumed to be specified, and the dust surface density is 0.01 
of the gas surface density; (iii) gas vertical structure is deter- 
mined from the hydrostatic equilibrium, while dust vertical 
structure is a consequence of turbulent stirring and grain set- 
tling. Also, calculating the chemical evolution we keep the 
dust properties fixed both for pristine and evolved dust cases. 
Below we describe each part of the model in detail. 

2.1. Radiative transfer 

The radiation in a protoplanetary disk plays a two-fold role. 
First, it is a main energy earner that redistributes energy com- 
ing from the stellar irradiation and viscous dissipation, and 
thus defines the overall disk structure. Second, it determines 
rates of photoreactions and thus shapes the disk chemical 
structure and observational appearance. These two aspects 
pose different requirements to the radiation transfer model. 
The radiation field as a contributor to the disk energy balance 
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IR wavelengths scattering can be considered negligible com- 
pared to absorption/emission. The p parameter describing 
anisotropy of scattering (p > 1/2 and p < 1/2 denote forward 
and backward scattering, respectively) can be introduced in 
our RT model in such a way that dust extinction efficiency a 
is substituted by the combination 2{\-p)a. In the limit of 
predominantly forward scattering grains, the role of UV dust 
heating in deep disk layers renders less significant than in the 
case of the isotropic scattering used in our study. That is, 
our current approach tends to slightly overestimate the role of 
scattering and thus overall dust heating in disk upper layers. 

Equations (l}-© are closed with the energy balance equa- 
tion 



4tt / K v {z)B v {T d {z))dv = A-k / nAz)JAz)dv + T m {z). (3) 



Figure 1. Overall computational scheme for ANDES. 



should be known in a wide range of wavelengths, from FUV 
(radiation from the accretion region and non-thermal radia- 
tion from the central star) to visual (thermal stellar radiation) 
to the infrared and submillimeter wavelengths (thermal disk 
radiation). This requirement makes multi-dimensional RT ap- 
proaches with high spectral resolution too slow for iterative 
disk modeling. 

On the other hand, in a narrow range of UV wavelengths 
(from 912 A to, say, 4000 A) such a good spectral resolu- 
tion is important for accurate calculation of the photochemical 
rates, as the dependence of photoreaction cross-sections on A 
is complicated. Protoplane tary disks usually have high optical 
depths at A < 100 ^m (e.g jBeckwith & Sargent! 19911) . which 
calls for using suitable methods to solve the radiation transfer 
(RT) problem for optically thick media. As our primary focus 
is on the chemical modeling in disks with evolved dust, we 
developed such a method with a particularly good wavelength 
resolution in the UV part of the spectrum. 

2.1.1. Main equations 

It is easy to show that in the cases of Schwarzschild- 
Schuster and Eddington approximations the RT equation for 
a plane-parallel ID medium can be written using the mean 
intensity J v : 



d 



1 dJ u (z) 



Xu(z) dz 



■■J v {z)-S v (z), 



(1) 



where xA cm -1 ] is the extinction coefficient, S v is the source 
function, and q = 1/4 and q = 1/3 for the Schwarzschild- 
Schuster and the Eddington approximation, respectively. 

If we consider only dust continuum absorption, thermal 
emission, and coherent isotropic scattering, the source func- 
tion is 

„ , , K u (z)BAT d (z)) + aJz)Uz) 

Sy(z)= — — — — . (2) 

KAz) + a u (z) 

Here ^[cm -1 ] is the absorption coefficient, a v [cm -1 ] is the 
scattering coefficient (Xu - K v +a v ) and B v is the Planck func- 
tion. 

In the 1+1D approach the anisotropic scattering by dust 
grains can also be taken into account. It is important for 
UV photons interacting with small dust grains, whereas at 



Here r nr (z) [erg cm" 3 s _1 ] accounts for non-radiative heat- 
ing/cooling mechanisms (gas-grain interaction, see Equa- 
tion (1A141 0. 

Equations (HJ-© represent the complete system for J„(z) 
and Ti(z). We s olve this system with the analogue of the 
Feautrier method (Mihalas] |1978l) . Specifically, we introduce 
a wavelength and coordinate grid where J„(z) is defined, and 
linearize the Planck function, B u , with respect to Equa- 
tion (HJ is approximated by a set of finite difference equations 
for each z-grid point, while Equation (01 is represented by a 
finite sum. As a result, we get a system of linear equations 
for J Vi (zk) that can be written using a hypermatrix formal- 
ism. This hypermat rix system is so lved with the tridiagonal 
Thomas algorithm dPress et al.ll992l) . After the new values of 
Jvi(zic) and 7d(Zfc) are obtained we refine linearization for the 
Planck function, update the system, and repeat iterations until 
convergence is achieved. 

The stellar and diffuse interstellar radiation fields can be 
treated as boundary conditions to the abov e system of equa- 
tions. We use an approach developed by Dull emond et al.l 
(120021) and consider stellar and interstellar fields as non- 
radiative additional source terms in Equation (0. This ap- 
proach takes into account the shielding of the star by the in- 
ner parts of the disk. For that one needs to know the frac- 
tion of stellar radiation intercepted by the disk at each ra- 
dius. We compute the corresponding grazing angle as an an- 
gle between dust density isoline at = 5 • 10~ 24 gem" 3 and 
the direction toward the star. For the stellar spectrum, we 
use a 4000 K blackbody for A > 4000 A . For shorter wave- 
lengths, we use the inte rstellar radiation field (tDrainelll978t 
Draine & Bertoldi 1996) with an ex tension to longer wave- 
length (van Dishoeck & Black 1982), where we have scaled 
the intensity so that it is continuous at the transition wave- 
length of 4000 A . Such a normalization leads to typical val- 
ues of stellar UV intensit y at disk atm osphere being equal 
to ~ 500 "Draine units" (R ollig et al.l 120071) at a radius of 
100 AU. 



2.1.2. Dust opacities and size distributions 

As a result of dust evolution modeling we get dust size dis- 
tribution functions /(a,/?,z)[cm" 4 ] being the fraction of grains 
with sizes within (a,a + da) interval. To compute dust opac- 
ities one should know efficiency factors for dust absorption 
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gabs and scattering gsc a : 

a max 



K v = I ira Q 3bs (a,u)f(a)da. 



a v = ira Q sca (a,v)f(a)da. 



(4) 



(5) 



(2 a bs an d Gsca ar e computed from th e Mie theory for as- 
trosilicate grains (lLaor & D raine 1993), but any other opacity 
model can be easily adapted. 

2.2. Gas thermal balance 

The kinetic gas temperature T g is obtained by solving the 
thermal balance equation: 



E 

k 



Tk(T g , T d , Pi ) - MT g , T d , Pi ,nf) = , 



(6) 



where T and A are gas heating and cooling rates in 
ergs -1 cm" 3 . They depend on absolute abundances of main 
heating/cooling species p, and their level populations ru, 
which in turn depend on the gas temperature. Therefore, 
the problem is solved iteratively at each grid point, starting 
from the disk atmosphere boundary toward the m idplane for 
any g iven radius, by means of the Brent method (Pre ss et all 
[19921 

Stellar FUV radiation is the main gas heating source in pro- 
toplanetary disks, leading to a PDR-like structure of the upper 
disk regions. There, gas is mainly heated via the photoelectric 
(PE) effect on dust grains and PAHs. In addition, collisional 
de-excitation of H2 pumped by FUV photons, photodissocia- 
tion of H2, and carbon photoionization are important heating 
sources in specific disk regions. Gas heating by exothermic 
chemical reactions plays only a minor role, with the largest 
contribution coming from H2 recombination on grains. In 
the optically thick, dense disk interiors, the dominant heating 
sources are the cosmic ray ionization of H and H2, and viscous 
heating due to dissipation of accretion energy. Gas mainly 
cools via non-LTE atomic and molecular line emission, colli- 
sions with grains, and, at high temperatures, by emitting Lya 
and metastable line emission. The details of heating and cool- 
ing processes can be found in Appendix lAl 

2.3. Chemistry 

An important ingredient of the thermal balance calculations 
is chemistry. While often a fast, simplified equilibrium ap- 
proach is adopted, time-dependent chemical modeling may 
be more appropriate for calculations of abundances of major 
molecular coola nts. We adopted the s ame gas-grain chem- 
ical model as in IVasvunin et al.l (1201 ll) . The reactions and 
reaction rates are ba sed on the RATE' 06 chemical rate file 
(IWoodall et al.l l2007). For all photochemical reaction rates, 
we use the local mean intensity (as a function of v) computed 
with the RT model. To compute photoreac tion rates, the disso- 
ciation and ionization cross-sections from lvan Dishoeck et al.l 
(2006) and the AMOP database^] are utilized. If cross-sections 
are not available for a certain reaction, we retain the standard 
Xoexp(-7Ay) formalism, with a 7 value taken from RATE'06 
ratefile, xo estimated at the upper disk boundary, and Ay 

^http: //amop. space . swri . edu/| 



computed as ln(x/xo)- The same values are used to esti- 
mate photodesorption rates. Thus, the calculation of photo- 
processes takes into account the detailed shape of the incident 
UV spectrum of the central star and its penetration through 
the dis k. Self-shielding for H? d issociation is computed us- 
ing the iDraine & Bertoldi] (119961) formalism, with the modi- 
fied Ay values used to account for dust attenuation. The self- 
and mutual shielding for C O photodissociation are computed 
using recent tabular data of Visser et al. (2009|). 
The unattenuated cosmic ray (CR) ionization rate is as- 



sumed to be 1.3 x 10 



The surface reactions are taken 



from lGarrod & Herbst (2006) and assumed to proceed with- 
out hydrogen tunneling via potential wells of the surface sites 
and reaction barriers. Thus only thermal hopping is a source 
of surface species mobility. A diffusion-t o-desorption en- 
ergy r atio of 0.77 is adopted for all species (Ruf fle & Herbstl 
2000). Under these conditions, stochastic effects in grain sur- 
face chemistry are negligible, and classical rate equations may 
be safely used (IVasvunin et alJl2009b iGarrod et al.ll2009l) . As 
the initial abunda nces, we u t ilize a set of "low metals" neutral 
abundances from lLee et al.l ([1998), where most of refractory 
elements are assumed to be locked in dust grains. 

As the density and temperatur e distributions, com puted 
here, are similar to those used in Vasy unin et al.l (1201 ll) . we 
decided to use the same vertical distributions of X-ray ioniza- 
tion rates regarding them as reference values. In the chemical 
module they are simply added up to CR ionization rates. For 
the purpose of chemical evolution, we assume that dust is rep- 
resented by grains with a single size which is computed from 
the loc al grain size distribution as described in IVasvunin et al.1 
(1201 ll) . With this chemical model, a disk chemical structure is 
computed using dust properties and physical conditions from 
the previous iteration. We assume that the grain properties do 
not change over the computational time span of 2 Myr. 

2.4. Vertical gas distribution 

Given that the gas temperature T g (R,z) and the mean molec- 
ular weight n(R,z) are known, the vertical gas density distri- 
bution can be found by integrating the hydrostatic equilibrium 
equation: 

9P(R,z) . GM,z 

T. =~p(R, Z)~ -372, (7) 

dz (R 2 + z 2 )' 

coupled to the equation of state of the ideal gas: 



kT a (R,z) 
P(R,z)= ^ ' \ p{R,z), 



(8) 



In this study we assume that the radial profile of the surface 
density is given by the known function £(/?). 

2.5. Dust evolution 

The evolution of the dust si ze distribution is cal culated us- 
ing the model presented in Birn stiel et al.l (12010b . In this 
work, we consider neither the viscous evolution of the gas 
disk nor the radial evolution of the dust surface density. 

The grain evolution begins with grains sticking at low col- 
lision velocities. Disruptive collisions at higher impact veloc- 
ities cause erosion or fragmentation, which poses an obstacle 
towards further grain growth and replenishes the population 
of small grains. Typical threshold collision velocities for the 
onset of fragmen tation are found to be about 1 m s" 1 for sil- 
icate dust grains (Blum & Wurm 2008). Icy particles are ex- 
pected to fragment at higher velocities due to the increased 
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surfac e binding energies dWada et al.ll2008t iGundlach et all 
120111) . We therefore use a threshold velocity for fragment- 
ing collisions of 10 and 30 m s" 1 in our dust models. Grain 
collisions are driven by re lative velocities due t o Brownian 
motion, turbulent motion (lOrmel & Cuzzil 120071) . radial and 
azimuthal drift as well as vertical turbulent settling. 

In order to make the calculation of the dust evolution feasi- 
ble, we consider a radially constant, vertically integrated dust- 
to-gas ratio and an azimuthally symmetric disk. The vertical 
settling of dust is ta ken into account by using a vertically in - 
tegrated kernel (see Braueret al. 2008; Birn stiel et al.l t2010). 
The integration implicitly assumes that the vertical distribu- 
tion of each dust species follows a Gaussian distribution with 
a size-dependent scale height. This is a good approximation 
for the regions close to the disk midplane where coagulation 
is most effective. However, for modeling of the chemical evo- 
lution the detailed vertical distribution of each dust species 
needs to be known. We th erefore use a vertical mix ing/settling 
equilibrium distribution (Dull emond & Dominikll2004h . tak- 
ing a vertical structure of the previously calculated dust sur- 
face densities. 

3. RESULTS 

3.1. Disk physical structure for evolved and pristine dust 

cases 

As an initi a l app roximation, we adopt a disk from 
iVasvunin et al.1 (1201 11) with mass M d i s k = 0.055M Q and gas 
surface density profile S(^) close to a power-law with in- 
dex p = -0.85 and S(1AU) = 34 g/cm 2 . The dust surface 
density is equal to 1% of the gas surface density. We as- 
sume the following parameters for a central star: a mass 
M* = 0.7M Q , a radius R+ = 2.64K Q and an effective tempera- 
ture = 4000 K. This system closely resembles the DM Tau 
disk. As UV-e xcess we u se the standard interstellar diffuse 
radiation field (Draine 1978) and its ex tension to longer wave- 
leng ths (v an Dishoeck & Blac k! 119821) as describ ed in Sec- 
tion |2Tl] ("JD" case from Akimk inet all (1201 ll) ). To show 
the influence of dust evolution on the disk thermal and density 
structure we present results for two cases: the pristine dust 
with uniform distribution and grain size of 0.1 /im (Model A) 
and the dust distribution and sizes as obtained with the dust 
evolution model after 2 Myr of evolution (Model Ev). The 
maximum grain size, attained in the midplane in Model Ev, 
varies from 4 • 10" 3 cm at 550 AU to 0.02 cm at 10 AU. The 
minimum grain size is always 3 • 10~ 7 cm. 

In Figure [2] the dust temperature distribution is shown for 
the both cases. The disk model with evolved dust is hotter 
by about 70 K in the inner disk atmosphere (R < 200 AU) 
and by ~ 10-20 K in the outer atmosphere (R > 200 AU) 
compared to the disk model with the pristine dust, whereas 
the dust midplane temperatures are quite similar in the both 
cases. Higher dust temperatures at the disk surface in the 
evolved dust model are explained by a steeper slope of dust 
opacities in the mid-IR, where such dust predominantly emits. 
Since both disk models have similar dust midplane tempera- 
tures and due to transparency of the bulk disk matter to the 
far-IR/millimeter emission, the emergent disk spectral energy 
distributions (SED) are similar. The difference in emergent 
spectra between Model A and Model Ev becomes apparent 
mostly at mid-IR wavelengths, where dust continuum emis- 
sion from the inner disk parts peaks. 

Gas temperature distributions in the disk models with the 
evolved and the pristine dust are shown in Figure [3] and can 
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Figure 2. Dust thermal structure for the disk model with the evolved (left 
panel) and pristine well-mixed (right panel) dust. 

be compared with the dust temperatures in Figure |2] The ex- 
tent of the gas-dust thermal coupling zone (where T& = T g ) in 
Model Ev is slightly smaller than in Model A, primarily due 
to sedimentation. As the midplane dust temperatures for the 
two models do not significantly differ, the gas temperatures 
also inherit this behavior. On the other hand, the gas tem- 
perature distributions above the coupling zone are quite dif- 
ferent. In the both cases, the inner disk atmosphere is heated 
up to several thousand Kelvin by photoelectric heating, but it 
is > 1000 K cooler in Model Ev. This is due to the reduced 
abundance of grains in the Model Ev, where the main contri- 
bution to the photoelectric heating comes from PAHs. 

In contrast, in Model A grains dominate photoelectric heat- 
ing. Their intense heating in the upper atmosphere has to be 
balanced by Lya cooling, while in Model Ev remaining grains 
in the inner disk and the [O I] line cooling at larger distances 
(> 40 AU) can balance the photoelectric heating from PAHs. 
Radial extent of hot tenuous atmosphere is drastically differ- 
ent for the two disk models: it exceeds 100 AU in Model A, 
whereas in Model Ev gas is cooler than 1 000 K even at R = 60 
AU. Absence of grains in the disk atmosphere in Model Ev 
leads to an increase of the gas temperatures by about factor of 
2 at z/R between w 0.3 and 0.6 at R « 100-300 AU. Some- 
what smaller increase of the gas temperatures in the upper lay- 
ers is also present in more distant regions of the disk. Rates 
of main heating and cooling processes are shown in Figure 2] 

Dust density distributions for Models Ev and A are com- 
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Figure 3. Gas thermal structure for the disk model with the evolved (left 
panel) and pristine well-mixed (right panel) dust. 

pared in Figure [5] The dust densities differ by several orders 
of magnitude, mainly due to sedimentation. This process also 
leads to dramatic changes in the dust-to-gas ratio. While in 
Model A this value is constant and equal to 10~ 2 , in Model Ev 
we encounter the whole range of values, from 10 -1 to 10~ 8 
(see Figure [SJ. However, dust-to-gas ratios below 10~ 4 lead 
to unstable solutions and poor convergence of the code, there- 
fore in the calculations we assume that the minimal value of 
dust-to-gas ratio is 10~ 4 . In Figure[7]the dust cross-section per 
hydrogen atom is presented for Model Ev. In case of Model A 
this value is equal to 5.9 • 10~ 22 cm 2 /H everywhere in the disk. 

3.2. Chemical structure 

One of the main goals of our study is to probe potential 
changes in the disk chemical structure that may arise due to 
various processes related to the dust evolution. In the section 
we present a detailed comparison of molecular abundances 
in the disk models with pristine and evolved dust for radii 
of 10 AU, 100 AU, and 5 50 AU. These a r e the same radii 
that have been analyzed by IVasyunin et al.1 (1201 ll) . We con- 
sider only those species that have mean abundances exceed- 
ing 10~ 12 at least in one of the two considered models. The 
mean abundance is computed as a ratio of the species column 
density to the column density of hydrogen nuclei. Remember 
that in all three cases we ignore and do not show the vertical 
structure at the height where the mass density of dust grains 



drops below the adopted limit of 5 ■ 10~ 24 g cm" 3 in the evolved 
model, and the medium can be considered as purely gaseous. 
In the case of the well-mixed dust disk model, this value cor- 
responds to the hydrogen number density of 2 • 10 2 cm" 3 , and 
is even lower for the dust evolving disk model. 

The key disk properties at the selected radii are shown in 
Figure [8] Up to a certain height the gas density is almost the 
same in both models. Above this height the vertical density 
profile flattens in Model Ev, because gas temperature either 
stops increasing or decreases with z, and density stays nearly 
constant to keep the disk hydrostatically stable. 

The main reason for that is the disk transparency. The dash- 
dotted lines in Figure [8] (b, d, f) show the water photodissoci- 
ation rates that are used here as a descriptive characteristics of 
the radiation field strength. Obviously, in Model Ev the UV 
radiation penetrates deeper into the disk. Dust is warmer in 
this model than in Model A almost everywhere in the disk. In 
Model Ev gas is also significantly hotter that in the model with 
pristine dust in the more illuminated region that extends ap- 
proximately from 1.5 AU to 3 AU at R = 10 AU, from 30 AU 
to 50 AU at R = 100 AU, and from 100 AU to 600 AU at 
R = 550 AU. 

We characterize the dust evolution using the total dust sur- 
face area per unit volume that is shown in Figure [8] (b, d, 
f) (solid lines). It is smaller in the evolved model as both 
grain growth and sedimentation reduce the total surface of 
dust grains. While in the midplane this reduction is mostly 
caused by the growth of dust grains and is quite moderate, 
from an order of magnitude at 10 AU to less than a factor of 
2 at 550 AU, in the upper disk, where sedimentation is impor- 
tant, the total dust surface area in Model A is greater by orders 
of magnitude than in Model Ev. However, this difference may 
not necessarily be important for chemistry as it is mostly con- 
fined to the illuminated disk regions where dust mantles are 
evaporated anyway by the UV photons. 

3.2.1. Outer disk 

We start the description of the disk chemical structure from 
the outer disk, where only minor changes in the disk physi- 
cal parameters are caused by the grain evolution. Because of 
ineffective grain growth the total grain surface area per unit 
volume is nearly the same in both models, except for the out- 
ermost disk atmosphere (Figure [8] f). Second, the dust tem- 
perature is quite low in the disk midplane, so surface reactions 
with heavy reactants should be mostly suppressed there. 

After 2 Myr of evolution we end up with 91 gas-phase 
species and 81 surface species that have mean abundances 
greater than 10~ 12 either in Model A or in Model Ev. In 
most cases grain evolution increases column densities for gas- 
phase components. Among the 91 gas-phase species only 13 
have column densities that are smaller in Model Ev than in 
Model A. The reason is quite straightforward. As grains grow 
and settle down toward the midplane, the so-called warm 
molecular layer moves down as well, to a denser disk region. 
Even if relative abundances do not change significantly in the 
process, column densities grow due to higher volume den- 
sities. Vertical abundance distributions for some species are 
shown in Figure [9] 

In the upper row of Figure |9]we present vertical abundance 
profiles for H2, H3, and CO. The key difference between 
Model A and Model Ev is that in the former the warm molecu- 
lar layer is located below the H2 dissociation boundary, while 
in the latter a portion of the molecular layer is located above 
this boundary, where free H atoms are abundant. This mutual 
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Figure 4. Main gas heating and cooling processes at selected radii. To avoid confusion, the legend is split into two parts. Same line styles apply on all plots. To 
illustrate the thermal coupling zone extent the cooling and heating functions are not plotted in regions with 7^ = Tg . 
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Figure 5. Distribution of the dust density in Model Ev (left panel) and Model 
A (right panel). 

disposition is not impossible as the molecular layer and the H2 
photodissociation front are not directly related to each other. 
Grain absorption of the FUV photons responsible for the H2 
dissociation is less significant in the model with evolved dust, 
and the H2 dissociation front is located deeper in Model Ev. 
However, the overall transparency of evolved dust in the en- 
tire UV range is smaller than for FUV photons only. Because 
of that the molecular layer that is controlled by desorption 
from dust grains and dissociation of trace molecules is located 



Figure 6. Dust-to-gas ratio for Model Ev. 

somewhat higher. This specific result, of course, depends on 
the adopted description of dust opacities and photoreaction 
rates. 

Simple atomic and diatomic components dominate the list 
of the gas-phase species, whose column densities are en- 
hanced in Model Ev. The largest column density increase 
at 550 AU is found for Si0 2 (Figure |9] d), N 2 0, and water 
(Figure [9] e). In all cases it is related to higher abundance 
of a molecule in the molecular layer. Note that almost all 
physical characteristics of the medium are nearly the same in 
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Figure 7. Dust cross-section per hydrogen atom for Model Ev. The corre- 
spondent value for Model A is 5.9 • 10 -22 cm 2 /H. 

the molecular layers of the two disk models, except for the 
gas density and the X-ray ionization rate. Higher density in 
the molecular layer of Model Ev accelerates two-body pro- 
cesses and shifts equilibrium abundances of many molecules 
to higher values. 

The relative location of the molecular layer and the H- 
H2 boundary, mentioned above, also plays a role in defining 
molecular column densities, especially, for species that are 
produced in reactions with atomic hydrogen, like water. In 
Model Ev a significant portion of the molecular layer is lo- 
cated above this transition, where abundant H atoms are avail- 
able. This speeds up the gas-phase water synthesis in H + OH 
reaction as well as surface synthesis of water molecules that 
are immediately released into the gas-phase due to photodes- 
orption. This explains a huge water spike located at height 
of about 300 AU. In Model A water is mostly produced in 
surface reactions that are less effective because of low H gas- 
phase abundance and lower rate of the UV photodesorption. 
Also, water molecules are more rapidly destroyed in reactions 
with ions such as HCO + that are abundant in the molecular 
layer due to higher X-ray ionization rate. The sharp drop-off 
in water abundance in both models coincides with the carbon 
ionization front. Above the front, the main destruction routes 
for water molecules are the reaction with C + and photodisso- 
ciation. 

The situation is somewhat different for complex hydrocar- 
bons, in particular, for long carbon chains. Their abundances 
in Model A are significantly enhanced in the molecular layer 
in comparison with Model Ev. This is again related to a more 
elevated position of the molecular layer in Model A. Because 
of that, it is less protected not only from the UV irradiation 
but from X-rays as well. Accordingly, ionized helium is more 
abundant in the molecular layer of Model A than in the molec- 
ular layer of Model Ev. Abundant C-bearing molecules, like 
CO, are destroyed by He + more efficiently in Model A in the 
disk upper region. Then, C + is consumed to produce sim- 
ple CH+ species that stick to grains and produce long carbon 
chains by surface processes. The dust temperature of the order 
of 30 K is high enough to drive desorption of these molecules 
into the gas-phase. 

This effect should not be overestimated. Even though the 
total column densities of carbon chains are greater in Model A 
than in Model Ev, their absolute values are low, with the 
mean abundance exceeding 10~ 12 only for C2H, C4, C4H, 
C5, C5H, and CeH The effect is most pronounced for C5H 



(Figure [9] f), with the ratio of column densities in Model A 
and Model Ev of 15. For the observed C2H molecule the 
higher relative abundance in Model A (related to more ef- 
fective He + chemistry) is nearly compensated by the higher 
absolute abundance in Model Ev (related to deeper location 
of the molecular layer), so its column densities are nearly the 
same in both models. 

The behavior of surface species is different in the two disk 
models. While column densities of gas-phase species are in- 
creased by grain evolution, column densities of many surface 
species decrease. There are 81 abundant surface species at 
550 AU, and only 30 of them have greater column densities 
in Model Ev. Also the difference of column densities of sur- 
face species is quite modest in the two models. Only for ten 
of 81 species column densities differ by more than a factor of 
3. Dominant surface carbon compounds (in terms of column 
densities) in both models are carbon monoxide and methane. 
Because of low dust temperature, S-CO2 production is sup- 
pressed, and this molecule in neither model reaches the high 
abundance seen at smaller radii (see below). 

Surface species that have greater column densities in 
Model Ev are mostly silicon and phosphorus compounds, 
which are not involved in surface chemistry (relevant reac- 
tions are not included in our chemical network). Their abun- 
dances are enhanced in the 'main' molecular layer as are 
abundances of their gas-phase counterparts (Figure |9] g). 

Abundances of some surface carbon chains are enhanced in 
Model A by about an order of magnitude due to more intense 
X-ray ionization than in Model Ev (see above). Also, species 
like S-C2O (Figure|9] h) and S-C2N involved in surface chem- 
istry have greater column densities in Model A because their 
midplane abundances are higher in this model due to greater 
available surface area for their synthesis. Carbon chains not 
involved in surface chemistry in our chemical model, like s- 
C4N (Figure|9] i), mirror evolution of their gas-phase counter- 
parts and have higher abundances in the upper carbon chain 
layer described above. 

3.2.2. Intermediate disk 

As we move closer to the star, at distances of about 50- 
100 AU, the fingerprint of dust evolution becomes more pro- 
nounced. While the mass density of dust is greater in the mid- 
plane of Model Ev due to sedimentation, the total surface area 
is still 2.5 times less than in Model A. In the upper disk the 
difference in the surface area reaches a factor of 70. It is in- 
teresting to note that the uppermost disk atmosphere is actu- 
ally colder in Model Ev than in Model A, despite being more 
transparent (Figure [8] c). This is because dust is not only the 
main source of opacity but also an important heating agent 
(due to photo-effect). As dust is depleted in the upper disk, 
the equilibrium temperature shifts to lower values, dictated 
by the PAH heating. 

At R = 100 AU, among 78 gas-phase species, having mean 
abundances higher than 10~ 12 at least in one of the two 
models, most species (72) have higher column densities in 
Model Ev, as at R = 550 AU, but the list of these species is 
somewhat different. Some examples of vertical abundance 
profiles for gas-phase species at R = 100 AU are shown in the 
top and middle row of Figure [TOl 

The main features of the chemical structure are the same as 
at 550 AU. In Model Ev the molecular layer, as marked by the 
CO distribution (Figure [10] c), is located above the H2 pho- 
todissociation front (Figure[l0] a). In Model A the situation is 
the opposite. Also, in Model A ions, like H3 (Figure[l0] b) are 
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Figure 8. Vertical distributions of selected disk parameters at three radii selected for a chemical analysis as indicated in titles. In all plots red lines correspond to 
a model with pristine dust, while blue lines correspond to a model with evolved dust. Shown in the left column are gas densities (dashed lines), dust temperatures 
(dotted lines), and gas temperatures (solid lines). Plots in the right column show dust surface area per unit volume (solid lines) and water photodissociation rates 
(dot-dashed lines). 



somewhat more abundant in the molecular layer which further 
decreases abundances of neutral unsaturated molecules. 

The largest difference between the two models is again ob- 
served for S1O2 that has a column density 3.6 • 10 8 cm" 2 in 
Model A and 1.8 • 10 12 cm -2 in Model Ev. Grain evolution 
causes water column density to increase from 7.1 • 10 13 cm" 2 
to 5.2 • 10 16 cm" 2 . This is again related to the different ar- 
rangement of the molecular layer and H-H2 transition in 
Model A and Model Ev (Figure [lOl a). The upper boundary 
of the water layer is defined by the location of the C ionization 



front. 

Many complex gas-phase hydrocarbons, like formaldehyde 
(Figure [TO] d) and cyanoacetylene (Figure [10] e), are also 
affected. Among more or less abundant molecules the only 
exception to this rule is methane (Figure [10] f), with col- 
umn density being about 3 times larger in Model A than in 
Model Ev. This difference is related to the surface chemistry 
as we will explain below. 

The chemical evolution of surface species is complicated 
as it is affected by at least two competing factors related to 
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Figure 9. Vertical abundance distributions of selected species at R = 550 AU. 



grain growth. Less grain surface is available to the mantle 
formation in Model Ev, but because of higher grain tempera- 
ture surface species are more mobile than in Model A, which 
intensifies the surface recombination. The interplay between 
these processes causes various responses of surface species to 
the grain evolution. 

Thirty five of 72 abundant mantle components have higher 
column densities in Model Ev. Only for 24 species the ra- 
tio of column densities in the two models exceeds a factor 
of 3. Higher abundances in Model Ev are mostly typical for 
complex surface molecules with large desorption energies that 
have accreted from the gas phase and are not involved in sur- 
face chemistry. These species mirror the behavior of their 
gas-phase counterparts. Two striking examples of greater col- 
umn densities in Model A are presented by s-CO and S-CH4 
(Figure [TO] g and h). Surface methane is underabundant by 
about 3 orders of magnitude in Model Ev, while surface CO 
is underabundant by more than 4 orders of magnitude in this 
model. As surface species are mostly concentrated in the mid- 
plane, to explain this difference we need to consider chemistry 
in this disk region. 



Surface methane chemistry is quite limited in the adopted 
network. Methane is synthesized in a sequence of hydrogen 
addition reactions, while the only effective route of its re- 
moval from the mantle is desorption. So, its abundance is 
regulated by the balance between hydrogen addition and des- 
orption. As the desorption energy of methane is not very large 
(1300K), desorption wins in this competition, and the steady- 
state S-CH4 abundance in Model Ev shifts toward lower val- 
ues. This does not work for other saturated molecules (like 
water and ammonia), as they have much higher desorption en- 
ergies, so their midplane abundances in both models are just 
(nearly) equal to the abundances of the corresponding atoms. 

Carbon monoxide is different from methane in the sense 
that it is not a 'dead end' of some surface chemistry route, but 
rather an intermediary on the route to S-CO2 formation and 
synthesis of complex organic molecules like CH3OH. The re- 
action between s-OH and s-CO, leading to S-CO2, has a 80 K 
barrier. This implies very strong dependence on the dust tem- 
perature at critical < 20-40 K. Because of this depen- 
dence, S-CO2 formation is much more efficient in Model Ev. 
While S-CO2 is the most abundant carbon compound in the 
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Figure 10. Vertical abundance distributions of selected species at R = 100 AU. 



disk midplane in both models, in Model A its abundance is 
about two times lower than in Model Ev (Figure[l0] i). Thus, 
in Model A, carbon atoms that are not bound in S-CO2 are 
available for other surface processes and are almost equally 
distributed between surface methane, surface CO, and some 
other species, with S-CH4 having almost the same abundance 
as S-CO2. In Model Ev, S-CO2 synthesis proceeds faster, and 
this component becomes not only the dominant C carrier, but 
also the dominant oxygen carrier, leaving almost no place for 
either surface or gas-phase CO. 

The described trends are mostly preserved at R = 50 AU. 
At this radius, dust evolution also causes column densities 
of abundant gas-phase species to increase (mostly due to en- 
hanced photodesorption). Water, carbon dioxide, formalde- 
hyde, and cyanoacetylene are among species mostly affected. 
Gas-phase methane column density is nearly the same in 
both models. Surface methane and CO ice are underabun- 
dant in Model Ev at R = 50 AU, but to a less degree than at 
i? = 100 AU. Dust temperatures in Model A and in Model Ev 
are nearly equal at this radius, so that in both cases surface 
S-CO2 synthesis is very effective, decreasing s-CO and S-CH4 



abundances and leveling differences between the two models. 
Abundances at 50 AU will be further discussed in the last sub- 
section of Section 3. 

3.2.3. Inner disk 

At R = 10 AU we have 75 abundant gas-phase molecules, 
and 62 of them share the common trend to be more abundant 
in the model with evolved dust. However, the magnitude of 
the difference in column densities as well as its origin are re- 
lated to other factors. The molecular layers both in Model A 
and in Model Ev are located above the H-H2 transition (Fig- 
ure QT| a). In both cases abundant H atoms are available 
both for surface and gas-phase reactions. Despite the fact, 
water column density in Model Ev exceeds that in Model A 
by more than 4 orders of magnitude (Figure [TT] g). This dif- 
ference is much greater than at other radii where we related 
it to the difference in atomic hydrogen abundance. At these 
warm temperatures of 50 - 200 K (see Figure [8] a), the for- 
mation of water is dominated by neutral-neutral reaction of O 
with H2 (with the barrier of 1 660 K), followed by the neutral- 
neutral reaction of OH with H2 (with the barrier of 3 163 K), 
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Figure 11. Vertical abundance distributions of selected species at R = 10 AU. 




h) 



see lWoodall etaD d2007h and lNajita et all d20TTl) . 

As H abundances in the molecular layers are nearly the 
same in both models, we need to find another explanation for 
the raise in water abundance in Model Ev. It is obviously 
related to the difference in physical parameters in the two 
molecular layers. Again, the molecular layer in Model Ev 
is shifted toward the midplane and, thus, resides in a denser 
disk region. Because of higher density in the molecular layer 
of Model Ev, surface water synthesis is more effective there, 
increasing its gas-phase abundance as well. Higher X-ray ion- 
ization rate in the molecular layer of Model A leads to higher 
ion abundances. In particular, it greatly enhances abundance 
of H + that is one of the major water destroyers. Another differ- 
ence is the UV radiation spectrum that favors carbon ioniza- 
tion in Model A. In the adopted photoionization model, car- 
bon atoms are ionized by the UV radiation with wavelengths 
shorter than 1 100 A. This radiation is absorbed less efficiently 
in Model A, and because of that the C/C + transition zone is 
further vertically expanded, so that the water layer falls in the 
region where C + abundance is still significant (Figure [TT] d). 
This also leads to rapid water destruction. 
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Different water abundances cause even greater differences 
in column densities of SO and SO2. In the case of SO2 the 
difference exceeds 9 orders of magnitude (Figure [TT] e). Sig- 
nificant growth of SO and SO2 abundances can be traced to 
the greater abundance of O2 in Model Ev. An SO2 molecule 
is produced from SO, SO is produced in reaction S + OH, 
atomic sulfur is the product of SO + dissociative recombina- 
tion, and SO + is produced in the reaction between S + and O2. 
Abundance of molecular oxygen in Model Ev is greater by 
almost 4 orders of magnitude than in Model A (Figure [TT] h), 
which is also related to different H + abundances, as the H + 
+ O2 reaction is one of the major O2 destruction pathways. 
Thus, we conclude that at R = 10 AU chemical differences be- 
tween Model Ev and Model A arise because grain evolution 
shifts the molecular layer in the region of the disk that is more 
protected from X-rays and FUV radiation. 

Among species, that have their column densities decreased 
by grain growth, are HCN (Figure [TT] f) and HNC. They 
are mostly concentrated in the midplane, and their midplane 
abundances in Model A exceed those in Model Ev by an or- 
der of magnitude. Analysis of chemical processes indicates 
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that this difference is related to surface processes, that is, 
higher gas-phase HCN abundance in Model A simply reflects 
more effective surface synthesis of the molecule because the 
available surface area is greater in this model (Figure QT] i). 
Then, HCN desorbs into the gas-phase and gets protonated by 
reactions with HCO + or H3. Dissociative recombination of 
HCNH + produces either HCN or HNC, so the overabundance 
of HCN in Model A is partially transferred into the overabun- 
dance of HNC. 

As for surface species, at this radius there are 64 abundant 
surface components, mainly heavy molecules. Nearly half of 
them are more abundant in Model Ev, but the increase in col- 
umn densities is not significant for most molecules. Two ex- 
treme examples of greater column densities in Model A are 
represented by s-HCN and s-HNC, for the reasons described 
above. 

3.3. Model with more efficient dust growth 

To check the sensitivity of our results to some details of 
the adopted grain physics, we considered additional models 
for dust evolution. In this subsection we present a detailed 
description of Model Evx with a threshold velocity for frag- 
menting collisions increased from 10 to 30 m s" 1 , which leads 
to more significant grain growth in the dense regions. 

Models A, Ev, and Evx can be viewed as successive stages 
of the grain evolution process. So, in Model Evx we may ex- 
pect to see a continuation of the same trends as were noted 
above for Model Ev. In Figure Q~2] we show the main disk 
structural properties at R = 50 AU in the three models. Ob- 
viously, more advanced grain evolution causes the disk to be- 
come more transparent. As a result, hot atmosphere becomes 
more extended, and dust becomes warmer, with the midplane 
grain temperature raising from 28 K in Model A to 33 K in 
Model Evx. As we will see below, this relatively small dif- 
ference has a noticeable effect on the disk chemical structure. 
Dust surface, available for chemical reactions, is an order of 
magnitude smaller in Model Evx than in Model A. 

Vertical profiles of some species for Models A, Ev, and Evx 
are shown in FigureQ~3] As in previous cases, we start from H2 
(Figure Q~3] a) and notice that the H2 photodissociation front 
sinks even deeper, so that hydrogen is almost fully atomic 
above ~ 6 AU. Due to warmer dust, gas-phase abundances of 
some molecules with low desorption energy are increased in 
the midplane of Model Evx (like in Model Ev at R = 10 AU). 
One can see the progressive growth of CO midplane abun- 
dance from Model A to Model Evx in Figure Q2](c). Similar 
to CO, gaseous N2 appears in the disk midplane in Model Evx. 
Protonation of such abundant molecules lowers the H3 abun- 
dance in the Evx model midplane (FigureQ~3] b), which affects 
abundances of some other ions, like H3O" 1 ". 

A typical example of the molecular abundance evolution is 
shown in Figure Qj] (d). A peak of water abundance shifts 
toward the midplane and grows higher. Due to increasing 
overall gas density and more intense photodesorption, gas- 
phase water column density increases up to 4.7 • 10 17 cm" 2 in 
Model Evx. The upper boundary of water layer is defined by 
the C + ionization front (Figure [T3] e). 

A significant growth is observed for N2H 1 " column density. 
It increases from 4.8 ■ 10 9 cm" 2 to 1 .5 • 10 10 cm" 2 in Model Ev 
and 5.1 • 10 11 cm" 2 in Model Evx (vertical abundance profile 
is shown in Figure |T3] (g)). This is related to increased ther- 
mal desorption of the N2 ice and lower abundances of surface 
species that are mostly synthesized on grains, like methane 



(Figure Qj] f) or ammonia (Figure Q~3] h). In the latter case, 
some nitrogen atoms in the midplane are free to be incorpo- 
rated into N2 molecules (Figure Q~3] i) and further into N2H + 
molecules. Species, significantly affected by the advanced 
grain growth, also include other gas-phase molecules, related 
to surface chemistry. Column densities are increased by more 
than an order of magnitude in Model Evx relative to Model Ev 
for H2O2, CH4, CO2, and some others. 

We have also considered the effects of dust radial mix- 
ing. The radial mixi ng is modeled as d i ffusio n, using the 
Schmidt number from lYoudin & Li th wield (12007). i.e. Z)d ust = 
f gas/ (1 + St 2 ). The dust diffusivit y is taken to be the dust vis- 
cosity (D gas = is g - ds ) which is the Shaku ra & Sunvaevl (11973b 
viscosity for the given alpha value. We found that radial mix- 
ing does not change the disk physical and chemical structure 
significantly and leads to the gas/dust temperature increase by 
several K at intermediate radii. 

4. DISCUSSION 

4.1. Comparison to \Vasvunin et al.\ ( 120771) 

While many aspects of the presented mo del are derived 
from the model used by Vasy unin et al.1 (1201 ll) . the new treat- 
ment of the disk structure results in parameters that are too 
different to allow a direct comparison of the 'old' and 'new' 
results. While density profiles are nearly the same in both 
studies, there are two key differences in the disk dust temper- 
ature and in the disk radiation field (Figure [14). First, the im- 
proved radiation transfer model makes dust in the 'new' disk 
midplane significantly warmer than dust in the 'old' disk mid- 
plane, at least, in the remote parts of the disk (R > 100 AU). 
Second, because of scattering the 'new' disk is less transpar- 
ent to dissociating far-UV radiation than the 'old' one. These 
two differences are related to the radiation tra nsfer treatment, 
so we may expect that basic inferences of IVasyunin et al.1 
(1201 ll) on the disk chemical structure should be retained in 
the new results, if they are mostly related to the dust evolu- 
tion. 

This is indeed the c ase, with a few except ions. First, the 
general conclusion of Vasvu nTn et al.l (1201 lb that dust evo- 
lution increases gas-phase column densities of most species 
is entirely confirmed in the present study. Second, almost 
all species, designat ed as sensitive to grain evolution in 
IVasyunin etail (12011b , like CO, C0 2 , H 2 0, C 2 H, retain this 
status in the present stud>0. 

In Table [l we show column den sities for species listed in 
Table 2 from lVasyunin et al] (1201 ll) . along with the newly cal- 
culated column densities. Few remarks are needed. Some 
species, like methanol, cyanopolyynes or formic acid, are sig- 
nificantly less abundant in the new model. This is due to 
generally less effective surface chemistry in a warmer disk 
midplane, where depletion of CO and other similar volatile 
ices is less severe and where a desorption rate of hydrogen 
atoms from dust surfaces is higher. The chances for them 
to be observed are, thus, slim (within the framework of our 
modeling approach). For some species from this list the 'sen- 
sitivity region' (the region where the two models differ most) 
shifts or extends to other radii (typically, from ten to hundred 
AU). These are HCNH + (derived mainly from HCN), NH 3 , 
and OH. Surface hydrogenation also plays an important role 
in the synthesis of these species. 

2 A species is marked as sensitive if its abundances in models with pristine 
and evolved dust differ by more than an order of magnitude at least at one of 
the considered radii. 
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Figure 12. Main disk structural properties at R = 50 AU. In all cases red lines correspond to Model A, blue lines correspond to Model Ev, and green lines 
correspond to Model Evx. In the left panel gas temperature (solid lines) and dust surface area per unit volume (dashed lines) are shown. Plotted in the right panel 
are dust temperature (solid lines) and gas number density (dashed lines). 
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Figure 13. Vertical abundance distributions of selected species at R = 50 AU. Blue dashed line corresponds to the Model Evx. 



Column densities of three species, H2CS, HC5N, and 
HCO + , while still enhanced by the dust growth, differ by less 
than an order of magnitude in the new calculation, so they 
do not conform to our sensitivit y criterion. Th i oform alde- 
hyde that has been mentioned in Vasv unin et al.l (1201 lb as a 
molecule most sensitive to dust growth and HCO + are now 
significantly more abundant in the midplane of Model A due 
to higher dust temperature and less severe depletion of their 
parental species, CO and (H)CS. This shortens the break be- 



tween the pristine and evolved dust models. 

Abundance of CH3CH3 is also significantly enhanced in the 
midpla ne at R = 10 AU, relative to results of IVasvunin et al.l 
(2011), and is nearly the same both in Model A and in 
Model Ev. As our model has a warmer inner midplane, sur- 
face radicals out of which CH3CH3 is formed become more 
mobile and reactive. A molecular layer no more dominates in 
its column density, so the molecule loses its sensitivity to the 
dust growth in the inner disk. In the outer disk the situation 
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z,AU z, AU z,AU 

Figure 14. Top row: disk structure in the present study and in Vasyunin et al. 12011), their Model GS. As in Figure[8] dashed lines show the gas density, and 
dotted lines show d ust temperature. Note that dust and gas temperature are equal in Model GS. Bottom row: water dissociation rates at the selected radii in the 
present study and in Vasyunin et al. ( 201 1). Different colors correspond to the same models as in the top row. 



is more complicated. There, CH3CH3 is still sensitive to dust 
growth, but the sign o f the sensitivity is different. While in 
Vasyunin et al. (2011) its column density was greater in the 
model with pristine dust, now CH3CH3 shares the common 
behavior and is enhanced in the molecular layer of Model Ev 
due to higher density there. 

Another molecule that shows the 'reversed' sensitivity is 
HCN. As we mentioned above, its abundance in the midplane 
is higher in Model A because of more effective surface syn- 
thesis. It also exceeds HCN column density in Model A5 from 
IVasvunin et al.l d201 ll) due to somewhat higher dust tempera- 
ture, that also intensifies HCN ice production (as surface pro- 
duction of CN is faster in the warmer ANDES model). At 
larger radii, HCN behavior is similar in both studies. These 
findings demonstrate the importance of the correct treatment 
of the radiation transport and also imply that the stellar and 
interstellar radiation fields need to be discretized as good as 
possible. 

Table Q] contains only a few representative species. To have 
a broader perspective, we perform a general comparison relat- 
ing the column density rati os in the model s with pristine and 
evolved dust computed in IVasvunin et al.l (1201 ll) and in the 
present study, for all species at 10, 100, and 550 AU. Results 
of comparison for R = 100 AU are shown in Figure [15] Only 
gas-phase species with mean abundances greater than 10~ 10 
are shown. Most species are concentrated around a red line 
that corresponds to equal old and new ratios. This indicates 
that m ost column densities re spond similarly to dust growth 
both in IVasvunin et al.l (1201 ll) and in the present study. Also, 
most species reside in the upper right quadrant, showing that 
in both studies dust evolution, as a rule, increases molecular 
column densities. 

Carbon dioxide is most sensitive to dust growth and is, thus, 
located in the upper right corner. This is not surprising as its 
production mostly occurs on grain surfaces via slightly en- 
dothermic reactions of CO and OH. A quite high water sen- 
sitivity was obtained in IVasvunin et al.l (1201 ll) . and now it be- 




10 -1 1Q o 1Q 1 1Q 2 1Q 3 

New evolved/pristine ratio 

Figure 15. R elation between molecular dust growth sensitivity in the present 
study and in Vasvunin et al. 12011). Plotted are ratios of column densities 
at R = 100 AU in the models with evolved and pristine dust. The red line 
corresponds to equal sensitivity. Species in upper right and lower left quad- 
rants show the same kind of sensitivity in both studies. Species within green 
lines have their column densities decreased by the dust growth in the previ- 
ous study, while in the present study dust evolution increases their column 
densities. A black square indicates difference in column densities less than 
an order of magnitude, which we interpret as an absence of strong sensitivity. 



comes even higher. Similar to water and carbon dioxide, HNO 
was very sensitive to dust growth in our old computation be- 
cause its main production route is surface synthesis. In the 
new computation this route is less effective due to warmer 
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Table 1 



Species sensitive to grain evolution. Observed column densities are compiled from Pietu et al. (2007), Dutrev et al. (2007), Chapillon et al. 12012), Bergin et al. 

<20TO|) . and Henning et alj )2010I)T 



Species 



CO 

co 2 

CH 3 OH 

H 2 

H 2 S 

C 2 H 

C 2 H 2 

CH3CH3 

H 2 CS 

HCN 

HC 3 N 

HC 5 N 

HC7N 

HCNH+ 

HCOOH 

OCN 

OCS 

NH 3 

HCO+ 

OH 



10 AU 



Pristine 



Evolved 



Column densities, cm 
100 AU 
Pristine Evolved 



550 AU 



Pristine 



Evolved 



2.0(17)6.5(17) 
6.8(10) 6.9(13) 
1.1(03)5.8(06) 
2.5(13) 1.3(12) 
1.7(05) 2.9(06) 
6.3(10) 1.8(11) 
6.2(10)3.1(12) 
3.0(09) 1.1(14) 
8.0(05) 6.4(09) 
2.6(12) 8.7(14) 
6.8(08) 6.8(07) 
4.9(08) 6.7(06) 
4.7(06) 5.8(02) 
2.0(10)2.1(11) 
1.4(11)6.6(06) 
3.2(07) 3.0(06) 
1.6(07) 1.9(06) 
1.7(11)3.7(10) 
1.1(11)2.7(12) 
6.1(13) 1.9(13) 



1.1(18) 1.8(18) 
9.5(13) 7.3(14) 
3.5(08) 2.6(08) 
1.1(16) 2.3(16) 
1.3(10)7.8(06) 
2.0(12) 1.0(13) 
4.6(12) 1.2(13) 
1.5(12) 6.9(14) 
3.5(11) 1.8(10) 
5.0(13) 4.5(13) 
1.7(12)7.5(11) 
1.5(12) 1.9(07) 
1.2(12) 2.1(06) 
2.7(11)2.1(11) 
6.3(13) 4.5(08) 
4.3(09) 5.5(07) 
9.3(10) 9.6(10) 
3.2(13)7.3(12) 
5.2(12) 8.1(12) 
1.9(14) 2.1(14) 



2.0(17) 
2.4(12) 
1.9(08) 
3.6(14) 
2.0(09) 
2.6(12) 
4.5(12) 
1.8(10) 
2.7(11) 
6.9(12) 
4.1(11) 
2.1(11) 
6.1(10) 
5.7(10) 
1.1(12) 
2.3(10) 
1.6(10) 
1.1(13) 
2.6(12) 
2.7(13) 



1.1(17) 
1.0(11) 
3.8(05) 
7.1(13) 
1.2(06) 
3.0(11) 
3.9(11) 
1.6(08) 
5.3(06) 
2.5(11) 
5.4(08) 
9.9(05) 
1.6(04) 
2.1(09) 
2.5(07) 
4.5(09) 
3.9(08) 
9.7(11) 
4.6(12) 
1.0(13) 



9.4(17) 
8.5(14) 
1.0(10) 
7.0(15) 
1.1(11) 
2.2(12) 
3.0(12) 
3.3(08) 
8.3(11) 
4.2(13) 
2.8(11) 
1.5(11) 
5.3(10) 
1.3(11) 
1.3(13) 
1.5(12) 
1.5(10) 
3.5(13) 
2.2(12) 
1.2(14) 



4.3(17) 
6.1(13) 
4.6(07) 
5.2(16) 
5.1(07) 
2.9(12) 
6.1(12) 
4.1(11) 
1.5(07) 
1.5(13) 
2.8(11) 
1.0(06) 
8.5(04) 
4.6(10) 
3.2(08) 
4.2(09) 
9.5(10) 
1.8(13) 
1.8(12) 
9.8(13) 



1.7(17) 1.3(16) 
8.2(13) 1.5(13) 
4.6(09) 1.2(07) 
1.4(14) 8.2(13) 
2.9(10) 6.3(09) 
6.7(12)5.7(12) 
4.4(12)3.0(12) 
1.3(07)5.1(05) 
2.7(10) 2.2(07) 
9.0(12)2.0(12) 
1.1(12)8.9(10) 
1.8(11) 1.6(07) 
9.2(10) 2.2(06) 
6.6(10)2.5(10) 
2.6(11) 1.3(09) 
1.3(11) 1.6(11) 
2.7(08) 1.4(08) 
8.5(12) 4.5(12) 
2.2(12)7.9(11) 
1.9(13)9.6(12) 



2.9(17) 
9.8(14) 
1.4(10) 
2.0(15) 
3.6(10) 
5.0(12) 
1.6(12) 
2.5(06) 
1.8(11) 
2.1(13) 
4.3(10) 
8.8(10) 
2.8(10) 
6.8(10) 
2.1(12) 
4.2(12) 
4.5(08) 
1.2(13) 
1.9(12) 
1-5(14) 



1.7(17) 
7.0(14) 
2.0(08) 
6.4(16) 
5.0(11) 
5.1(12) 
8.3(12) 
8.4(07) 
2.2(07) 
2.5(13) 
6.6(11) 
2.4(07) 
9.8(06) 
8.6(10) 
1.3(09) 
4.3(11) 
3.9(09) 
2.3(13) 
1.3(12) 
2.2(14) 



Observed column 
densities, cm -2 
DMTau 



1.4(17) 

< 3.0(13) 
2.8(13) 

6.5(12) 



6.5(12) 



dust, so HNO is mostly produced in the gas phase. This 
makes it less susceptible to changes in dust properties. It must 
be kept in mind that warmer dust has dual effect on surface 
chemistry. On one hand, a larger temperature implies more 
rapid hopping and larger reaction rates. On the other hand, 
more volatile reactants evaporate faster from warmer grains, 
thus, quenching the formation of some species. 

An opposite example is represented by formaldehyde. This 
species was barely sensitive to dust growth in the old com- 
putation, with column density being slightly smaller in the 
model with evolved dust. In the present study, H2CO column 
density is significantly greater in the model with evolved dust. 
This behavior is related to details of the UV penetration. In 
the old models, where only the UV absorption has been taken 
into account, the UV field intensity falls off quite slowly as we 
go deeper into the disk. Because of that abundance maxima 
in the molecular layer for molecules that are mostly suscep- 
tible to photodesorption and photodissociation are extended 
and shallow. Thus, their column densities are less sensitive 
to dust evolution. Detailed treatment of the radiation transfer 
in the new model predicts a sharper transition from the illu- 
minated atmosphere to the dark interior. The molecular layer 
becomes significantly narrower and is, thus, much more sen- 
sitive to the extent of dust growth and sedimentation. 

The overall conclusion from the presented comparison is 
the following. We confirm that dust evolution changes col- 
umn densities of many molecules (see Tableland Figure [TBI. 
Most species that have been listed as e specially sensitive to 
dust evolution in Vasvu nin et alj (1201 ll) retain this status in 
the present study. However, column densities of some species 
turn out to depend on the details of the radiation transfer treat- 
ment, and this dependence will become even stronger when 
we will proceed from column densities to line intensities. AN- 
DES makes all the necessary preparatory work for that, pro- 
viding us with both abundances and gas temperatures. 

However, there is another aspect, apart from the radiation 
transfer, that may affect our conclusions. This aspect is re- 



lated to possible evolutionary changes. As in ANDES we use 
time-dependent chemistry, we can provisionally estimate its 
importance. 

4.2. Disk structure with time-dependent chemistry 

In order to demonstrate the effect of time-dependent chem- 
istry on the disk chemical structure we perform model calcu- 
lations with abundances of major molecular coolants at 10 4 , 
10 5 , and 2 x 10 6 years. We assume that the disk chemically 
evolves from mostly neutral atomic gas, with molecular hy- 
drogen and a low fraction of atomic hydrogen (10~ 3 to the to- 
tal number of hydrogen nuclei). We do not consider the time 
evolution of dust grain distribution in order to focus on effects 
of chemical evolution, and utilize vertical dust distribution af- 
ter 2 Myr. Results are presented in Figure Q~6] showing the 
relative abundances of H2, H, CO, C, and C + as a function 
of height at the distinct epochs for disk radii of 10, 100 and 
550 AU. 

As can be clearly seen, the location of the H2/H transition 
shifts toward the midplane with time for all the considered 
radii due to slow photodissociation of molecular hydrogen, 
self-shielded from the strong FUV stellar radiation. H2 cannot 
be quickly re-formed in this region in Model Ev due to overall 
lack of grains, providing catalytic surface for H + H reaction. 
Consequently, between 10 4 and 210 6 years, at radii of 10, 
100, and 550 AU, the PDR zone shifts from 1.6 to « 1.4 AU, 
from 31 to 22 AU, and from 425 till 290 AU, respectively. 
This effect is more pronounced in the outer disk, where den- 
sities and density gradient are lower. 

An interesting feature of the chemical structure in Model 
Ev is the presence of a 'dip' in H2 vertical distribution at the 
final time moment at all radii. This region with depression 
in H2 concentration is caused by its slow X-ray/UV destruc- 
tion, which cannot be compensated by the H2 surface pro- 
duction on a few remaining grains. However, just above this 
depression region dust-to-gas ratio locally increases, and so 
does the available surface for hydrogen recombination (per 
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unit gas volume). The reason for the elevated dust-to-gas ra- 
tio is a gas redistribution from the top of the coupling region 
to greater heights due to extra heating. 

In contrast, the evolution of ionized carbon reaches a chem- 
ical steady-state rapidly everywhere in the disk thanks to its 
fast ion-molecule chemistry pathways, so the C + concentra- 
tion is not time-sensitive (after 10 4 years), see Figure Q~6] The 
chemical evolution of C + is governed by a simple and limited 
gas-phase reaction network in the disk atmosphere, where it is 
an imp ortant coolant with a relative abundance of ~ 1 0~ 4 (see, 
e.g. ISemenov et"ai1 12004 ISemenov & Wiebd 12011b . Neu- 
tral atomic carbon shows little time evolution, if any, in the 
disk regions adjacent to the midplane at radii smaller than 
~ 100 AU. Due to relatively large densities in disk equatorial 
regions, neutral carbon is rapidly converted to CO and hydro- 
carbons. This is not true for lower-density outer disk regions, 
at R > 500 AU and z/R ~ 0.3-0.6, where the C abundance 
changes substantially with time. Since initially all elemental 
carbon is in the atomic form, in the outer disk, less dense and 
more transparent to the interstellar FUV radiation, conversion 
of C into CO and hydrocarbons takes more than 10 4 years (see 
Figure [16] bottom row). 

The gas-phase CO abundances follow the pattern of H2/H 
and C and do not reach a steady-state within 2 Myr every- 
where in the disk model with evolved dust. The grain growth 
increases the CO freeze-out timescale to > 1 Myr in the in- 
ner and intermediate radii (Figure [16] top and middle rows). 
In the midplane, where gas-phase CO abundance is low, this 
molecule is present as CO ice. The final distribution of the 
CO abundance at R < 100 AU shows an interesting feature: 
due to severe grain growth CO freeze-out is inefficient in the 
midplane, but still effective at disk heights of z/R ~ 0.2 AU. 
At even larger heights the CO molecular layer starts, so CO 
emission lines are excited both in the very cold and warm re- 
gions. Since 12 CO, 13 CO and C ls O isotopologue lines, hav- 
ing vastly different opacities, allow probing these two tem- 
perature zones, this should be visible with modern radio- 
interferometers. Intriguingly, evi dence for the presen ce of 
very cold CO gas was found by Dart ois et al.l (120031) . and, 
later, for other m olecules like HCO+, CC H, CN, and HCN 
(see discussion in Semenov & Wiebe 20llt). 

Enhanced amounts of H2 at 10 4 years in disk regions with 
high FUV radiation intensities lead to additional heating by 
the UV-excited H2. Since the H2/H boundary is moving down, 
the gas thermal structure of the disk responds accordingly and 
also shows strong variations of T g in a narrow disk layer, in 
particular, at R > 100 AU (see Figure \T7\. While the gas tem- 
perature varies from 250 K to w 200 K (25%) atR = 10 AU, at 
the outer disk region, R = 550 AU, the gas temperature differ- 
ence at various times is about 250 K (from 320 K to w 75 K, 
or a factor of 4), compare top and bottom panels of Figure [P71 
Naturally, it should also have a strong impact on chemical 
composition and appearance of the disk molecular layer, from 
which most of line emission emerges. More importantly, 
it demonstrates the importance of using the time-dependent 
chemistry for calculating abundances of key gaseous coolant 
instead of the commonly applied steady-state approach. 

5. CONCLUSIONS 

A multi-dimensional self-consistent model of protoplane- 
tary disks ANDES' is introduced and described. The purpose 
of ANDES is to provide researchers with a state-of-the-art, 
most up-to-date detailed thermo-chemical model of a proto- 



planetary disk that can be used to analyze high-quality IR and 
(sub-)millimeter observations of individual nearby disks. For 
the first time grain evolution and large-scale time-dependent 
molecular chemistry are included in modeling of physical 
structure of protoplanetary disks. 

The iterative ANDES code is based on a flexible modular 
structure that includes 1) a 1+1D continuum radiative transfer 
module to calculate dust temperature, 2) a module to calculate 
gas-grain chemical evolution, 3) a 1+1D module to calculate 
detailed gas energy balance, and 4) a 1+1D module that cal- 
culates dust grain evolution. The disk structure is computed 
iteratively, assuming fixed dust density structure after the first 
iteration. Typically it takes ~ 10 iterations to reach conver- 
gence at 1 % level of accuracy. 

The continuum radiative transfer module is based on the 
two-stream Feautrier method with a high-resolution fre- 
quency grid. We consider dust continuum absorption, ther- 
mal emission, and coherent isotropic scattering. The dust evo- 
lution is modeled by accounting for coagulation, fragmenta- 
tion, and gravitational sedimentation towards the disk mid- 
plane balanced by turbulent upward stirring. The chemical 
model is based on a gas-grain realization of the RATE' 06 net- 
work, and includes surface reactions and X-ray/UV processes. 
All modules have been thoroughly benchmarked with previ- 
ous studies, with overall good agreement and performance. 

We study the impact of dust evolution on dust tempera- 
ture, gas temperature, and chemical composition by compar- 
ing results of the disk models with evolved and pristine dust. 
We compute gas thermal structure corresponding to chemical 
abundances evolving from the initial abundances for 10 4 , 10 5 , 
and 2 • 10 6 years. We show that time-dependent chemistry 
is important for a proper description of gas thermal balance. 
The strongest impact on the gas temperature (up to 100 K) 
occurs in the outer, low-density region beyond 100 AU. This 
is mainly due to the shift of H2/H PDR transition deeper into 
the disk with time. 

In accordance with previous studies, it is found that the gas 
becomes hotter than the dust in elevated disk regions reaching 
1 000-10000 K in the inner atmosphere. However, the main 
heating source is different for the two dust models. In the 
disk with pristine dust it is photoelectric heating by grains. In 
the atmosphere of disk with evolved dust grains are strongly 
depleted, therefore photoelectric heating by PAHs becomes a 
dominant heating process. Thus a realistic, observationally- 
based estimates of absolute PAH abundances and sizes are 
required to calculate accurately gas temperature in the inner, 
~ 1 - 50 AU disk atmosphere accessible with Spitzer, Her- 
schel, and ALMA. 

The response of disk chemical structure to the dust growth 
and sedimentation is twofold. First, due to higher trans- 
parency a partly UV-shielded molecular layer is shifted closer 
to the dense midplane. Second, the presence of big grains in 
the disk midplane delays the freeze-out of volatile gas-phase 
species such as CO there, while in adjacent upper layers the 
depletion is still effective. Even though the dust evolution 
shifts the molecular layer of the water vapor closer toward 
the cooler, midplane disk region, it increases its overall con- 
centration. This aggravates the disagreement between the wa- 
ter vapor column densities predicted by modern astrochemi- 
cal models, which are higher than t hose observed with Her- 
schel in the disks around TW H ya dHog erheijde et al.ll20llh 
and DM Tau (Bergin et al. 2010) b y factors of at l east several 
(see also discussion in Semenov & Wiebd (1201 ll) ). Overall, 
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Figure 16. The left panels show H>/H transition at selected radii calculated for different epochs: 10 , 10 , and 2 ■ 10 yrs for the disk with evolved dust. The 
right panels show the CO/C/C + transition. 



molecular concentrations and thus column densities of many 
species are enhanced in the disk model with dust evolution, 
e.g., CO2, NH 2 CN, HNO, H 2 0, HCOOH, HCN, CO. 
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APPENDIX 

A. MAIN GAS HEATING AND COOLING PROCESSES 
A. 1 . Main heating processes 

PHOTOELECTRIC HEATING BY GRAINS We follow iKamp & van ZadelhofJd2001h and calculate the photoelectric heating rate by 
silicate grains as 

r PE = 2.5 x 10-\ u b Ie dustX , (Al) 
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where is the dust absorp tion coefficie nt at UV wavelengths, x is the strength of the UV radiation field in units of the 
Draine FUV interstellar field (Draine 1978), and e^ust is the photoelectric efficiency determined by the grain charge parameter 
x = yJT^x/n e (here n e is the electron number density). For e we adopt expressions from Kam p & van Zadelhof fl (12001b . The 
relative strength of the FUV, \, is defined as 

r 1 lOnrn , , , 
J91.2nm A ^A 

rllOnm . Draine . ' 
J91.2nm AM A " A 

The average dust opacity at UV wavelengths is determined by integration of frequency-dependent dust absorption cross-sections 
in the UV frequency range 

C = f [ f(aWQ, bs (a)dad Vl (A3) 



where f(a) is given by the dust evolution model. 

PHOTOELECTRIC HEATING BY PAHS Polycyclic aromatic hydrocarbons (PAHs) possess large cross-sections for UV photon ab- 
sorption and therefore can efficiently heat gas by photoelectric emission, even if their abundance is low. Heating by PAHs can 
be particularly important for disks with evolved dust, since PAHs are better mixed wi th gas than macroscopic d ust particles, and 
thus re main in the disk atmosphere while bigger grains settle towards the midplane (Dulle mond et al.l 2007a). Bake s & Tielensl 
( 119941) derived a simple analytical expression for their PE heating rate: 

= 10- 24 /pAH«HepAHX, (A4) 

where n H is the hydrogen nuclei number density, and e PAH = 0.0487/(1+4 x 10~ 3 x 073 ). The parameter / PAH is the depletion 
factor of the PAH a bundance relative to the diffuse ISM value, which is estimated to be ~ 10-20% of the total carbon budget 
dDraine & Lill2007l) . The details of the evolution of PAHs in protoplanetary disks are far from being fully understood, though it is 
clear that high-energy stellar radiation may play an enormous role in th eir destruction and chemical transformation ( Acke et ail 
2010; Siebenmorgen & Krugel 2010; Siebenmorgen & Hevmann 2012). Therefore, we do not consider PAHs in the simulations 
of dust evolution and treat /pah as a free parameter of the model. In the present paper we as sume /pah = 0-1 based on estimates 
from observations of PAH spectra in disks surrounding young T Tauri and Herbig Ae stars (Keller et al. 2008; Kamp l201ll) . A 
detailed study of the effects of PAHs heating on the structure of protoplanetary disks with evolved dust is beyond the scope of the 
present paper. 

COSMIC RAY HEATING Cosmic ray (CR) particles deposit energy mainly through ionization of H2 and H at the rate 
dBakes & Tielenslfl99l : 

Tcr = Ccr (5.5 x 10" 12 «(H) + 2.5 x 10- n »(H 2 )) , (A5) 
where Ccr s -1 is the attenuated CR ionization rate and n(X) denotes concentration of a species X. 

HEATING BY SURFACE H 2 FORMATION Formation of one H 2 molecule on the grain surface liberates 4.48 eV of energy, but the 
exact partitioning of this energy into H 2 vibration, rotation, translation and accommodation by a grain lattice remains uncertain. 
It is commonly assumed that this energy is equally redistributed between rotational, vibrational and translational movements. We 
assume that formation of one hydrogen molecule returns only 1 .5 eV (2.4 x 10~ 12 erg) to the gas (IB lack & Dalgarnoll976l) . Then, 
the corresponding heating rate is 

r H ,form = 2.4 X 10" 12 R H ,fonn«H, (A6) 

where/?H,f 0lm is the H 2 formation rate in s" 1 . The further details of calculation of chemical reactions rates are described in 
Sect. [Of 



PHOTODISSOCIATION OF H 2 We take into account only spontaneous radiative d issociation of H? : H? +hv - > H 2 * -> H + H + /ija 
Assuming that the average kinetic energy of dissociation products is 0.45 eV (Stephens & Dal garnolll973l) . the corresponding 
heating rate is 

r H2 dis = 6.4 x 10- 13 R H2 phdis«(H 2 ), (A7) 

where ^H 2 phdis is the p hotodissociation ra t e of H 2 . To calculate this rate, we take into account self-shielding of H 2 molecules as 
given by Eq.(37) from lDraine & Bertoldil (11996b . 

COLLISIONAL DE-EXCITATION OF H 2 In dense PDR regions col lisional de-excitation of FUV-pumped W 2 is the second 
most important heating mechanism (Stern berg & Dalgarnol Il99~5b . Here we adopt a simple two-level approximation of 
H 2 vibrational heating and coo ling from Rtillig et al. (2006), which nevertheless well reproduces the net heating rate computed 
bv lSternberg & Dalgarnol (1 19951) with 15 vibrational levels. The net vibrational heating is given by the following expression: 

r„ c , = r H .-A H2 , (A8) 

where Th* is the vibrational heating rate by collisional de-excitation and Ah 2 is the vibrational cooling rate. For details of the 
calculation of the heating and cooling rates we refer to the Appendix C in Rollig et al. (2006). 
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C photoionization — Ionization of atomic carbon releases electrons with kinetic energies of - 1 eV dBlacklll987l) . The corre- 
sponding heating rate can be approximated as: 

r c =1.6xl(T 12 Rcph«(C), (A9) 
where Rc p h is the photoionization rate of the C atoms. 
Viscous heating — The viscous heating rate is given by dFrank et al.lll992b : 

Tvis = ^kin^ ep , (A10) 

where the kinematic viscosity of the gas is parameterized as i>\± n = acjH g ( Sh akura & Su nvaev 1973), cj is the isothermal sound 
speed, H g is the gas pressure scale height, and Q,\^ v is the Keplerian velocity. 

A.2. Main cooling processes 

NLTE LINE COOLING The net line cooling rate for a given species is determined by the total amount of upwards and downwards 
radiative transitions. Level populations for each coolant are calculated from statistical equilibrium equations. Unlike FUV, the 
local FIR intensity that enters these equations depends on the temperature and level populations over the large part of the disk. 
This requires iterations over all vertica l grid points simultaneously . To simplify a calculation, we adopt an escape probability 
approach using the expression (B9) in Tiele ns & Hollenbachld 19851) . 

We perform the full non-LTE calculations, considering the major coolants for a typical PDR: fine structure lines of C, O, C + and 
rotational lines for the CO molecule. The data for energy levels, collision, emission and absorption coefficients for computation 
of the NLTE line cooling are taken from the LAMDA database (Schoier et al. 2005). The data include collision rate coefficients 
for collisions of H2 , H, e", He, and H + with O and C atoms, as well as collisions of H2 , H, and e" with C + , and H2 with CO. For 
minor coolants we use approximate formulas presented below. 

HIGH-TEMPERATURE COOLANTS The cooling by emission from metastable levels of neutral and ionic species becomes impor- 
tant at temp eratures exceeding several t housand Kelvin. We calculate the cooling rate from 'D- 3 P emission by O I (630 nm) 
according to Sternberg & Dalgarnol (119891) : 

A O i630 = 1.8 x 10- 24 n(O)« e exp(-22800/r g) , (All) 

where «o is the neutral oxygen concentration. 

The coolin g by electron impact ex citation of metastable levels of ionic species (e.g., Fe + , Fe" 1 " 1 ", Si + ) is calculated by approximate 
formula from lDalgarno & McCravl (11972b : 

A ion (T)=A i T-^ 2 e xp(-T i /T s ), (A12) 

where parameters A, and T, for each ion are given in Dalg arno & McCravl A 19721) . 

Another im portant cooling process at high temperatures is Lya emission. We adopt the cooling rate by Lya emission from 
lSpitzeri(fl978T ): 

A Lya = 7.3 x 10" 1 V«(H)exp(-118400/7 g ). (A13) 

H 2 LINE EMISSION Rotational line emission of the H2O molecule can contrib ute to cooling in dense disk regions. We include 
line c ooling rates o f H2O due to the excitation by H2 , using analytical fits from lNeufeld & Kaufmanl d!993l) for T > 100 K and 
from lNeufeld"etal] ( fl995l) for 10 K< T <100 K. 

THERMAL ACCOMMODATION Thermal accommodation is the energy exchange by inelastic collisions between dust and gas. In 
disk models with standard ISM-like dust it is a dominan t cooling process, with the exception of the upper, tenuous atmosphere 
and ou ter radii (R > 400 AU) (e.g.. iWoitke et al.ll2009t) . We utilize the corresponding cooling rate from Bur ke & Holle nbach 
(fl983l) : 

A d _ g = 4x 10- n Tr(a z )n d n H a T ^/T„(T g -T d ), (A14) 
where the thermal accommodation coefficient otj is set to 0.3 ( a typical value for silicates and carbon). 

B. RADIATIVE TRANSFER AND GAS THERMAL BALANCE BENCHMARKING 
The RT module of the ANDES code was checked for the following cases which allow an analytic or semi-analytic solution: 

• Optically thin case, zero non-radiative heating source (1A14I) . arbitrary incident radiation. In this case the mean intensity 
in the media is equal to the incident one and temperature may be derived from (01 and it is the same for every position in 
media. Tests show an equality of temperatures which are derived by our RT code and by independent numerical solution 
of the energy balance equation ||3}. 

• Arbitrary media (optically thick, non-gray opacities) with the non-diluted Planck incident radiation. In this case the media 
becomes isothermal with temperature of radiation. Test show an equality of temperatures and radiation field derived from 
analytical and numerical solutions. 
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• Optically thick media with gray opacities, arbitrary incident radiation and non-radiative dust heating function which is 
proportional to dust density and has no other z-dependence. In this case the system of Equations ([T])-© narrow down to 
the second order linear ODE and may be solved analytically. A comparison of the analytical and numerical solutions is 
presented in Figure[l8] The mean intensity parabolic profile is shown in the inset graph as well. 



^ 400 





RT solution 
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Figure 18. Computed and analytical temperatures and mean intensities for a test case. 




The surface of the disk is a photodissociation region (PDR) controlled by UV radiation from a star and interstellar radiation 
field. Therefore, w e perform benchmarking of the thermal balance in our code as proposed in the PDR code comparison study 
(Rollig et al. 2007). For benchmarking purp oses we use a red uced chemical network restricted to the most abundant elements 
(H, He, O, C, e~) and 31 species(Table 4 in Roll ig et alJl2007h . The reaction rates are taken from the UMIST99 database with 
some corrections from A. Sternberg. H? disso ciation rate is 5 x IO '^y/IOs" 1 . Cosmic ray H ionization rate is C = 5 x 10~ 17 s _1 . 
For more details of benchmark test we refer to Rollig et al. (2007h. 

All benchmark models assume plane-parallel, semi-infinite geometry of clouds of total constant hydrogen density of 10 3 and 
10 5 5 cm" 3 . The values of the standard far UV field were taken as x = 10 and 10 5 times the Draine (1978) field. There are two 
sets of benchmark models: four with fixed dust and gas temperatures of 30 and 50 K, respectively, and the other set of four 
models with the gas temperature resulting from thermal balance. The first set of models with fixed temperature aims at testing 
main ingredients of the thermal balance: solutions of chemistry and statistical equilibrium equations for level populations of main 
coolants, while the second set examines solution of thermal balance. Here we present results of benchmark tests for both kinds 
of models with density ng 1 = 10 5 5 cm" 3 and far UV field strength \ = 10 s (models F4 and V4 in lRollig etail (120071) ). 

The left panel of Figure [19] shows comparison of our calculations with post-benchmark results for the H/H2 transition zone 
typical for PDR environment. Right panel of Figure [l9]shows the C + /CO/C transition zone. 

Main heating and cooling rates included in benchmarking are shown in the left panel of Figure [20] Gas-grain cooling and 
[OI] 63/imline are the dominant cooling processes for A y < 0.5. CO l ines dominate cooling at high attenuated regions. Our line 
cooling rates show remarkable agreement with data from Rollig et al. (2007) for dominant cooling processes: [OI] 158/im, [OI] 
63, 145 /mi, [CI] 370, 610/imlines. Comparison of our model results for gas temperature in the slab with other PDR codes is 



Protoplanetary disk structure with grain evolution: 



THE ANDES MODEL 



25 




Ay Ay 

Figure 20. Left pa nel. Comparison of gas temperature calculated by our code (red solid curve) with post-benchmark results from different PDR codes from 
Rollig et al. 12007). Thick magenta dashed line marks average gas temperature from the benchmarking data. Right panel. Main heating and cooling processes 
included in our code for the benchmark model V4. 



shown in Figure [20] At small Ay the gas temperature is much higher than the dust temperature due to photoelectric heating and 
agrees well with other PDR codes. 



